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Non-Gaussianity of the primordial density perturbations provides an important measure to con¬ 
strain models of inflation. At cubic order the non-Gaussianity is captured by two parameters tnl 
and ^nl that determine the amplitude of the density perturbation trispectrum. Here we report mea¬ 
surements of the kurtosis power spectra of the cosmic microwave background (GMB) temperature as 
mapped by Planck by making use of correlations between square temperature-square temperature 
and cubic temperature-temperature anisotropies. In combination with noise simulations, we find 
the best joint estimates to be tnl = 0.3 zb 0.9 x 10^ and ^nl = —1.2 zb 2.8 x 10^. If tnl = 0, we find 
= -1.3 zb 1.8 X 10^ 

PACS numbers: 


Introduction .—Existing cosmological data from cosmic 
microwave background (CMB) and large-scale structure 
(LSS) are fully consistent with a simple cosmological 
model involving six basic parameters describing the en¬ 
ergy density components of the universe, age, and the 
amplitude and spectral index of initial perturbations. 
The perturbations depart from a scale-free power spec¬ 
trum and are Gaussian. These facts support inflation 
as the leading paradigm related to the origin of density 
perturbations mm- Under inflation a nearly exponen¬ 
tial expansion stretched space in the first moments of 
the early universe and promoted microscopic quantum 
fluctuations to perturbations on cosmological scales to¬ 
day m [5]. Moving beyond simple inflationary models 
with a single scalar field, models of inflation now involve 
multiple fields and exotic objects such as branes that have 
non-trivial interactions. Such inflationary models pro¬ 
duce a departure from Gaussianity in a model-dependent 
manner mi- The amplitude of non-Gaussianity there¬ 
fore is an important cosmological parameter that can dis¬ 
tinguish between the plethora of inflationary models [TO] . 

The first order non-Gaussian parameter, /nl, has been 
measured with increasing success using the bispectrum - 
the Fourier analog of the three-point correlation function 
of the CMB temperature. Such studies have found /nl to 
be consistent with zero EHni, with the strongest con¬ 
straint coming from Planck given by /nl = 2.7±5.8 [15]. 
The inflationary model expectation is that /nl ^ 1 and 
a constraint at such a low amplitude level may be fea¬ 
sible in the future with large scale structure data and 
with 21-cm intensity fluctuations. Alternatively, with 
the trispectrum or four point correlation function of CMB 
anisotropies [T6| , we can measure the second and third or¬ 
der non-Gaussian parameters tnl and ^nl- While these 
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higher order parameters generally lead to a trispectrum 
that has a lower signal-to-noise ratio than the bispec¬ 
trum, there may be models in which the situation is 
reversed with the trispectrum dominating over the bis¬ 
pectrum contribution. An example of such a model is 
an inhomogeneous end to thermal inflation discussed in 
Ref. HI]. 

A previous analysis using WMAP data out to < 600 
using the kurtosis power spectra involving two-to-two 
and three-to-one temperature correlations nHum, found 
-7.4 < ^nl/ 10^ < 8.2 and -0.6 < tnu/IO"^ < 3.3 
at the 95% confidence level (C.L.). Other measures 
of the WMAP trispectrum have been presented in [20]- 
1^ . While the Planck data have been used to constrain 
Tnl < 2800 at the 95% C.L. such a constraint ignored the 
signal associated with ^nl m- Using all of the Planck 
data, the expectation is that ^nl can be constrained with 
a 68% CL uncertainty of 6.7 x 10^ [21] with tnl = 0, 
while Tnl can be constrained down to 560 if ^^nl = 0 
[24] . Here we present an analysis of the Planck temper¬ 
ature anisotropy maps by making use of kurtosis power 
spectra to constrain tnl and ^nl jointly. 

Theory .— We begin the discussion with the tempera¬ 
ture trispectrum defined as [25] 


{CLI 

imi ^ l 2 m 2 ^^47714 ) 


f h h L \ 

^ ^ I mi m2 —My 

LM ^ ^ 

(i i 


where we have introduced the Wigner 3-j symbol. The 
angular trispectrum, (L), can be further expressed 
in terms of sums of the products of Wigner 3-j or 
6-j symbols times the so-called reduced trispectrum, 


tU:{l) US]. 

To derive the angular trispectrum given by Tj^l^ (L) we 
assume that the curvature perturbations ( of the universe 
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generated by inflation follow as: 


$(x) = $g(x)+/nl($g(x)-($g(x)))+c/nl$g(x) . (2) 


where the curvature perturbation ( and the initial grav¬ 
itational potential are related by ^ = (3/5) C and tnl = 


(6/nl/5)^ 

We refer the reader to Ref. [24] for intermediate steps in 
our derivation. Using the above form the full trispectrum 
can be reduced to two forms involving the two amplitudes 
Tnl (associated with 4>^(x) — (4>^(x)) term in above) and 
^nl coming from 4>^(x). 


Defining = 

1, 2 [26], where 
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FIG. 1: The estimator validation using WMAP simulations 
with Tnl = 3600. 
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we find that the reduced trispectrum is 

7;y=(L) = [rNL7;y-(')(L) + 5NL7/y-(')(L)]. (4) 

The two terms are 

= ^ nl (^) j rldrirldr2FL{ri,r2) 

aiAri)Mri)ai^ir2)l3uir2), (5) 

and 

= s-nl J r‘^drl3i^{r)l3^{r)[ai,{r)(3i^{r) 

+ah{r)^i,{r)]. (6) 

Here ctz(r) = {2/7r) J k‘^dkAf^{k)ji{kr) and 

^i{r) = {2/7r) f k‘^dkP{k)A^^{k)ji{kr). The pri¬ 

mordial curvature power spectrum is k^P{k)/{27r^) = 
{3/b)‘^As{k/ko)^^~^ with no “running” [27] • Here ko is 
the pivot scale set at 0.05Mpc“^. We use the public 
code [33] to compute ct/(r), A(r) and the temperature 
transfer function Aj^{k). 

In the Tnl part, we define the function Fl as 

FL{ri,r2) = ^ J k‘^dkP{k)jL{kri)jL{kr2). (7) 

Following the efficient algorithm in [28], we define ^ = 
r^jrx^x = kr\ and compress r\ and into one dimension 
such that 

J dxx^>-^jL{x)jL{tx), (8) 

Here A = (3/5)^(27r^//co)As/co~’^^ We validate that this 
fast algorithm gives the same results as Eq. 

The first part of the trispectrum associated with tnl 
approximates to y/Ci-^Ci^Ci^Ci^ at L < 100. 

This is due to the fact that the integrand peaks at r = 
and Cl = f r‘^drai{r)/3i{r) [29]. Here r* is the comoving 


distance at last scattering surface and C'^* = PL^r^^r^). 
For the comparison with the data, however, we perform 
an exact calculation defined in Eqs. mu The adaptive 
r-grid is used for the integration. 

The estimators of the connected trispectrum are con¬ 
structed in Refs. [liisn] and they are given by 


^i^’^^(^NL,^NL) 


and 


1 

tlTT 


E 

I 1 I 2 I 3 I 4 


1 

2L + 1 Ci^Ci^Ci^Ci^ 

(9) 


Uf’^l3-NL,5'NL) = 




dih 


2h 


1 , 2 -^ + 1 O 1 O 2 O 3 O 4 

I 1 L 2 I 3 X 


( 10 ) 

In Eqs. 9, 10, the reduced trispectrum Ti^i^ {L) is eval- 
uated at tnl = 1 and ^nl = 1- The estimators iF/ ’ ^ 
and are parametrized by these two parameters. 

The denotes the full trispectrum from data or 

simulation. 

In our analysis, /min E /15/25/35T E /max^ /min — 

2 and /max = 1000. The trispectrum computing 

time is proportional to 0(/max) at a single L. In 
order to make these calculations more efficient, we 
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use Monte Carlo integration for iF/’ y i.e., replacing 

E ^max \ ''^max \ ''^max \ ''^max Utt- TT / ]\T \ ^ 

Z^Z2=Zmin 2^l3=lmin ^ Samples Z^i- 

The vector l(=/i, / 2 ,/s,/ 4 ) is uniformly sampled from 

[/min? /max] aud V = (/max /min) • Eor K , WC re¬ 
strict the diagonal elements within 2 < L < 20 and val¬ 
idate that a bigger upper bound negligibly modifies the 
trispectrum. The Wigner 3-j symbols’ intrinsic selec¬ 
tion rule also helps reduce the computation time. With 
all these efficient algorithm, we can achieve a hour-level 
computation time, which is about three orders of magni¬ 
tude faster than the brute-force calculation. We show the 
theoretical predictions of these estimators for the case in 
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Fig.[T] for a fixed set of tnl and ^nl values for which 
non-Gaussian simulated maps are available. 

From simulated and real data, spherical harmonic 
coefficients and are computed by inverse 

spherical harmonic transformation (SHT). Then the two 
weighted maps are generated from definitions A{r, n) = 
'EimMr)aimYimi-a) , B{r,n) = Eim 
and aim = where the angular power spectrum 

Cl is inclusive of noise, is calculated by anafast of 

Healpix which removes monopole and dipole. To correct 
the masking effect, we scale the masked modes and 

^(^ata) py 2/ fsky to match the underlying temperature 
power spectrum. These masked modes are also beam- 
and pixel window-deconvolved. In the following text, we 
neglect “n” for brevity. 

From A and B maps, we construct C(ri,r 2 ) = 
A{ri)B{r 2 ). Then we make = ^^(ri,r 2 )Cz^(ri,r 2 ) 
and T)(ri,r 2 ) = G'(ri,r 2 )A(r 2 ). We can calculate four 
types of power spectra: 

= ^^A„^(rl,r2)BL(r2), (11) 

m 

jAb,aB(^i^^ 2) = _^^F((ri,r2) 

m 

[AB]Un)[AB]l^{r2)-, (12) 

Fy^’^r) = ^ ^[^BB],„(r)BL(r); (13) 

m 

and 

^ J2i^B]Ur)[BB]Ur). (14) 

m 

When all the power spectra are integrated along the 
line of sight, they become: 

jABA.B ^ y r^rfrir|dr2J;^®^’®(ri,r2); (15) 

LfBB.B ^ J ^2^^^ABB,B(^). 
jAB,AB ^ y r2drir|dr2J/^B’^B(ri,r2); (17) 

and 

LfB.BB ^ J 

The trispectrum estimators 

Kp) = (5yjAB,AB^2Ly’BB^ (19) 


and 
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are then constructed from the correlations associated 
with A and B maps that are either from data or sim¬ 
ulations. 

These estimators are applied to 143 GHz and 217 GHz 
temperature datasets, as well as the cross-correlation 
143 X 217 GHz. For the cross correlation, the estima¬ 
tors are 


k[^’^^(143 X 217) 


— tA(143)B(217),A(143)B(217) 

_l_ 2£^A(143)B(217),B(143)B(217) (2^) 


and 

ii'f4)(i43^217) = (5y^A(143)B(217)A(143),B(217) 


Simulation Validation: To validate our estimates of 
the connected trispectra, we make non-Gaussian GMB 
signal simulations. The non-Gaussian maps for WMAP 
are publicly available [34] so we simulate maps with 
^side = 512 and /max = 600, and all the WMAP experi¬ 
mental settings, consistent with 5-year observations, are 
adopted. For the signal part, aim = and we 

choose /nl = 50, i.e., tnl = 3600 given the expected rela¬ 
tion between /nl and tnl, independent of the exact value 
of ^fNL- Note that the non-Gaussian simulations we use 
assume ^nl = 0 and in a joint model fit to data we test 
this expectation. The WMAP 5-yr noises are then added 
in the signal simulations. The WMAP simulation is 
F(n) = biPiaimYim{n) + CTo/y77(n)nwhite(n). Here 
(Jo and A^(n) are provided by WMAP. The estimator of 
the connected trispectrum is Kl = l/4!(i^L — 

In Fig. 12 we show that the average connected parts from 
100 full-sky realizations are consistent with the theoreti¬ 
cal calculations. 

Data Analysis and Results: We use Planck 143 GHz 
and 217 GHz temperature maps for the present analy¬ 
sis. We use the foreground mask to remove the point 
sources and galactic emissions for both frequencies. The 
217 GHz map cleaned after the 70% foreground mask still 
contains visible emission around the galactic plane, so we 
use an extended mask to further cut the 217 GHz data 
around it. The resulting sky fractions for both maps 
become 73% and 58%. At 143 GHz, the map is con¬ 
volved with a 7' Gaussian beam and has 45/iK arcmin 
noise. At 217 GHz, it is 5' and 60/iK arcmin. Follow¬ 
ing Ref. m, point sources (PS) and cosmic infrared 
background (GIB) are also included in simulated data. 
The power spectra for these two sources are = 

27r/3000^ and = 27r/(/(/ -j- l))(//3000)^‘^, respec¬ 

tively. The foreground power at these frequencies are 
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FIG. 2: The raw trispectra calculated from Planck data and 
simulations for 143 x 143 GHz (top) and 143 x 217 GHz (bot¬ 
tom). In both plots Gaussian bias dominates the raw signal. 


with the parame- 

ters ^143x143 = 64/4X^,74^^3x217 = 43/4X^,7421^7x217 = 
57/4X2, AfigBx 

143 — 4/4X2, 

217 ~ 14/4X2, AC^B^' 

217 — 

54/iK^. In addition, a 10/iK arcmin white noise is added 
into the simulations. The data structure is expressed as 
= ^im^imkpiyimi'^) + ri{n) where n is a direc¬ 
tion on the sky, bi is the beam transfer function, pi is the 
pixel transfer function at rigide = 2048, and n(n) is the 
noise simulation. We use 100 signal and noise realiza¬ 
tions from the FFP6 simulation set of the Planck collab¬ 
oration [32] . We use the best-fit cosmological parameters 
from “Planck-fWP+highL” [27] • Specifically, = 

0.022069, nch‘^ = 0.12025, r = 0.0927, = 0.9582, 

^4^ = 2.21071 X 10“^ at pivot scale ko = 0.05Mpc“^, and 
Ho = 67.15km s“^Mpc“^ [27]. 

We calculate both trispectra and from 

Gaussian simulations and data for Planck. The Gaus¬ 
sian term in the trispectra is averaged from 100 

Planck simulations for frequency combinations 143 x 143 
GHz, 143 X 217 GHz and 217 x 217 GHz, and is removed 
from the raw signal, which is defined as the combination 



FIG. 3: The 68%, 95% and 99% confidence levels for different 
combinations are indicated by the transparency of the con¬ 
tours. The frequency combinations 143 x 143 GHz, 143 x 217 
GHz and 143 x 143 + 143 x 217 GHz are shown in blue, red 
and black colors. 

of the connected part and the disconnected part. All 
the trispectra are shown in Fig. It is seen that the 
disconnected components dominate the raw signal and 
our simulations can precisely recover these significant 
biases. Also, all the trispecta show consistent shapes. 
From 100 simulations, the full covariance matrix M is 
obtained for each frequency combination and the vec¬ 
tor Vh = Here b is index of trispec¬ 

trum band. We choose five bands for each spectrum: 
L=[2,152], [152,302], [302,452], [452,602], [602,800]. Here 
we use AL = 150 and Lent = 800. We want to both 
avoid systematic issues with the high L trispectra and 
get enough signal-to-noise, so we choose this conserva¬ 
tive cut here. 

We choose a binning function to maximize the sensi¬ 
tivity 

T> _ . C _ ^Leb 

here Sl = {2L ^ 1)Kl is the fiducial model with tnl = 
^NL = 1, A^L = (2A -b ^ -b 1)^^^ 

which is the connected trispectrum from the simulation 
or data. 

The likelihood function of the data is given as 
xb^NL,5NL) = 

bb' 

(24) 

where the two free parameters are tnl^^nl, b index of 
the band, and u the index of the frequency combination. 
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TABLE I: The constraints of tnlj^'nl with AL = 150 and 
Tcut = 800 from different frequency combinations. The 68% 
C.L. is given by = 2.3 except the last row. 

Freq.Combination tnl[x10^] ^nl[x10^] 

143 X 143 -0.6 ± 1.2 -1.9 ± 3.9 

143x217 1.9 ±1.5 -1.0 ±4.1 

143 X 143 ± 143 X 217 0.3 ± 0.9 -1.2 ±2.8 
143 X 143 ± 143 X 217 0 -1.3 ±1.8 


TABLE II: The constraints of tnlj^'nl with different AL and 
Lcut for the combination 143 x 143 ± 143 x 217. The 68% C.L. 
is given by = 2.3. 

143 X 143 ± 143 X 217 tnl[x10^] ^nl[x10^] 

[AL = 150, Lcut = 800] 0.3 ±0.9 -1.2 ±2.8 

[AL = 150, Lent = 850] 0.3 ± 0.9 0.3 ±1.5 

[AL = 150, Lcut = 900] 0.4 ±0.9 1.7 ±1.4 

[AL = 200, Lcut = 800] 0.6 ±0.9 -0.6 ± 3.0 


We draw O(IO^) samples for two parameters from 
Monte Carlo Markov chains with flat priors —10^ < 
^NL < 10^ and —10^ < ^nl < 10^. The 217 GHz map is 
still significantly contaminated by CIB although we use 
a very conservative cut which removes 40% of the sky, 
so we do not include 217 x 217 GHz into our parameter 
estimation. The constraints for tnl and ^nl are listed 
in Table [T[ In the last row of Table |Tj we show the 1- 
parameter constraint on ^nl with tnl = 0. For all the 
combinations, we find that tnl and ^nl are consistent 
with zero. We check the consistency between different 
frequency combinations in Fig. From Fig. it is seen 
that different bin sizes do not change the results. We also 
check the impact of effective L range on the parameters. 
From Figj^ we find that adding more L range can result 
in a higher value of ^nl and the interpretation is that 
the high L range is systematically contaminated by un¬ 
resolved point sources and non-Gaussian contribution of 
CIB beyond the foreground mask. All the results shown 
in Fig. are summarized in Table [Hj 

Summary: We present the first joint constraints 
on tnl^^nl using Planck kurtosis power spectra that 
trace square temperature-square temperature and cubic 
temperature-temperature power spectra. The Gaussian 
biases in these statistics are corrected for with simula¬ 
tions and we make use of non-Gaussian simulations to 
test our pipeline. We find the best joint estimate of the 
two parameters to be tnl = (0.3 ± 0.9) x 10^ and ^nl = 
(-1.2 ± 2.8) X 10^ If Tnl = 0, ^nl = (-1.3 ± 1.8) x 10^ 
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FIG. 4: The 68%, 95% and 99% confidence levels for the 
combination 143 x 143 ± 143 x 217 with different bin sizes 
(top) and Lcut (bottom) are indicated by the transparency of 
the contours. In the top, for AL =150, the contour is shown 
in black and green for AL = 200. For both cases, Lcut=800. 
In the bottom, Lcut = 800 is shown in black, Lcut = 850 in 
red, Lcut = 900 in blue. In these cases AL = 150. 
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